Conventionally, a known image restoration algorithm restores a blurred image, in which degradations have occurred in the image imaged by an imaging apparatus such as a digital camera, due to being out of focus, hand movement, aberration and the like.
As the image restoration algorithm, for example, a known method represents image degradation due to blurring on the captured image with a degradation function (point spread function (PSF)), and restores it to an image having no blurring based on the degradation function.
Image restoration algorithms using the degradation function are known including, for example, the Wiener filter, the general inverse filter, the projection filter and the like. Japanese Patent Laid-Open Publication No. 2004-205802 discloses an image restoration method using the Wiener filter, and Japanese Patent Laid-Open Publication No. 2004-288653 discloses an image restoration method using the general inverse filter.
Here, the degradation function used by the method of restoring the blurred image will be described.
When f(x,y) is an ideal image and g(x,y) is the blurred image, it is assumed that there is a relation as follows:
[Formula 1]g(x,y)=∫∫h(x,y,x′,y′)f(x′,y′)dx′dy′+v(x,y)  (1)
Here, h(x,y,x′,y′) represents the degradation function, and v(x,y) represents random noise in an output image.
Except for the parallel translation, if an image having a blurred point does not depend on a position of the point, the degradation function becomes h(x−x′,y−y′) and the expression (1) becomes as follows:
[Formula 2]g(x,y)=∫∫h(x−x′,y−y′)f(x′,y′)dx′dy′+v(x,y)  (2)
When there is no noise, the Fourier transformation is performed on both sides of the expression (2) and it becomes as follows according to the convolution theorem:G(u,v)=H(u,v)F(u,v)  (3)
G(u,v), F(u,v) and H(u,v) represent the Fourier transformations of g(x,y), f(x,y) and h(x,y), respectively. Here, H(u,v) is a transfer function of a system for transforming the ideal image f(x,y) into the blurred image g(x,y).
Now, in order to consider only the degradation due to the relative movement of the camera and the landscape, its degradation model will be described below. If the relative movement approximately equals that caused by the movement of an imaging area of an image sensor in a plane, the total of exposures at a point on the imaging area is obtained by integrating instantaneous exposures only when a shutter is opened. In this case, it is assumed that the time required for opening and closing the shutter can be ignored. When α(t) and β(t) are components of X and y directions of a displacement respectively, it becomes as follows:
[Formula 3]g(x,y)=∫−T/2T/2f(x−α(t),y−β(t))dt  (4)
In this case, T is an exposure time, which is assumed to be from −T/2 to T/2 for descriptive purposes. If the Fourier transformation is performed on both sides of the expression (4), it becomes as follows:
                    [                  Formula          ⁢                                          ⁢          4                ]                                                                                                                G                ⁡                                  (                                      u                    ,                    v                                    )                                            =                            ⁢                              ∫                                                      ⅆ                    x                                    ⁢                                      ∫                                                                  ⅆ                        y                                            ⁢                                                                                          ⁢                                              exp                        ⁡                                                  [                                                                                    -                              j                                                        ⁢                                                                                                                  ⁢                            2                            ⁢                                                                                                                  ⁢                                                          π                              ⁡                                                              (                                                                  ux                                  +                                  vy                                                                )                                                                                                              ]                                                                                                                                                                                                          ⁢                                                ∫                                                            -                      T                                        /                    2                                                        T                    /                    2                                                  ⁢                                                                  ⁢                                  ⅆ                                      tf                    ⁡                                          (                                                                        x                          -                                                      α                            ⁡                                                          (                              t                              )                                                                                                      ,                                                  y                          -                                                      β                            ⁡                                                          (                              t                              )                                                                                                                          )                                                                                                                                              =                            ⁢                                                ∫                                                            -                      T                                        /                    2                                                        T                    /                    2                                                  ⁢                                                                  ⁢                                                      ⅆ                    t                                    ⁢                                      ∫                                                                  ⅆ                        x                                            ⁢                                              ∫                                                  ⅆ                                                      yf                            (                                                                                          x                                -                                                                  α                                  ⁡                                                                      (                                    t                                    )                                                                                                                              ,                                                                                                                                                                                                                                                                                  ⁢                                  y                  -                                      β                    ⁡                                          (                      t                      )                                                                      )                            ⁢                                                          ⁢                              exp                ⁡                                  [                                                            -                      j                                        ⁢                                                                                  ⁢                    2                    ⁢                                                                                  ⁢                                          π                      ⁡                                              (                                                  ux                          +                          vy                                                )                                                                              ]                                                                                                    
If x−α(t)=ξ and y−β(t)=η, the above expression becomes as follows:
                                              ⁢                  [                      Formula            ⁢                                                  ⁢            5                    ]                                                                                                                        G                ⁡                                  (                                      u                    ,                    v                                    )                                            =                            ⁢                                                ∫                                                            -                      T                                        /                    2                                                        T                    /                    2                                                  ⁢                                                                  ⁢                                                      ⅆ                    t                                    ⁢                                      ∫                                          ∫                                                                        ⅆ                          ξ                                                ⁢                                                  ⅆ                          η                                                ⁢                                                                                                  ⁢                                                  f                          ⁡                                                      (                                                          ξ                              ,                              η                                                        )                                                                          ×                        exp                                                                                                                                                                                  ⁢                                                [                                                            -                      j                                        ⁢                                                                                  ⁢                    2                    ⁢                                                                                  ⁢                                          π                      ⁡                                              (                                                                              u                            ⁢                                                                                                                  ⁢                            ξ                                                    +                                                      v                            ⁢                                                                                                                  ⁢                            η                                                                          )                                                                              ]                                ⁢                                  exp                  ⁡                                      [                                                                  -                        j                                            ⁢                                                                                          ⁢                      2                      ⁢                                                                                          ⁢                                              π                        ⁡                                                  (                                                                                                                    α                                ⁡                                                                  (                                  t                                  )                                                                                            ⁢                              u                                                        +                                                                                          β                                ⁡                                                                  (                                  t                                  )                                                                                            ⁢                              v                                                                                )                                                                                      ]                                                                                                                          =                            ⁢                                                F                  ⁡                                      (                                          u                      ,                      v                                        )                                                  ⁢                                                      ∫                                                                  -                        T                                            /                      2                                                              T                      /                      2                                                        ⁢                                                            exp                      ⁡                                              [                                                                              -                            j                                                    ⁢                                                                                                          ⁢                          2                          ⁢                                                                                                          ⁢                                                      π                            ⁡                                                          (                                                                                                u                                  ⁢                                                                                                                                          ⁢                                                                      α                                    ⁡                                                                          (                                      t                                      )                                                                                                                                      +                                                                  v                                  ⁢                                                                                                                                          ⁢                                                                      β                                    ⁡                                                                          (                                      t                                      )                                                                                                                                                                  )                                                                                                      ]                                                              ⁢                                                                                  ⁢                                          ⅆ                      t                                                                                                                                              =                            ⁢                                                F                  ⁡                                      (                                          u                      ,                      v                                        )                                                  ⁢                                  H                  ⁡                                      (                                          u                      ,                      v                                        )                                                                                                          (        5        )            
From this expression (5), it can be seen that the degradation is modeled with the expression (3) or the expression (2) which is equal to the expression (3). Then, the degradation function H(u,v) is given by the following expression:
[Formula 6]H(u,v)=∫=T/2T/2exp[−2jπ(uα(t)+vβ(t))]dt  (6)
Therefore, the degradation function in the case of blurring with regard to the X-axis at an angle θ, at a certain speed V and only in the time T is given as follows:H(u,v)=sin πωT/πω  (7)
Here, ω is as follows:ω=(u−u0)V cos θ+(v−v0)V sin θ  (8)
u0 and v0 are central coordinates of the image. It should be noted that when ω is very small, the degradation function approximates to H(u,v)=T.
As the image restoration algorithm using the degradation function H obtained in this way, some methods have been proposed including a method of performing the Fourier transformation with regard to the blurred image and performing an operation in a frequency space, or a method of performing an iterative computation with a steepest descent method with regard to the image in a real space, and the like.
Here, the method of performing the Fourier transformation with regard to the blurred image and performing the operation in the frequency space will be described.
In other words, it is assumed that the blurred image g(x,y) and the original image f(x,y) are in accordance with the model of the above described expression (2). Then, if there is no noise, the Fourier transformations of g(x,y), f(x,y), and the blurred image H, h(x,y) satisfy the above described expression (3). Here, the expression (3) is transformed as follows:F(u,v)=G(u,v)/H(u,v)  (9)
From this expression (9), if H(u,v) is known, it is possible to restore f(x,y) by multiplying the Fourier transformation of the blurred image G(u,v) with 1/H(u,v) to perform the inverse Fourier transformation.
Also, as the other method, the method of performing the image restoration by using the iterative computation with the steepest descent method will be described by using FIG. 18.
FIG. 18 is a flowchart showing a processing procedure of an image restoration processing apparatus for performing the image restoration by using the iterative computation with the steepest descent method.
In FIG. 18, the image restoration processing apparatus sets the blurred image g(x,y) as the 0th restored image (that is, initial image) in a RAM (S10). Next, a parameter n indicating an iteration count is initialized to 0, and further a maximum iteration count nMAX is read from a ROM (S12), and a predetermined convergence parameter ε is read from the ROM (S14). Also, a threshold Thr is read from the ROM as an end determination parameter (S16). Next, if the iteration count n is less than a predetermined maximum iteration count (a determination result at step S18 is positive “Y”), after the iteration count n is incremented (S20), VJ (nabula) is calculated (S22), and a square of norm of VJ is calculated to be set as a parameter t (S24).
Here, J is an evaluation amount of the general inverse filter and given as follows with the blurred image g(x,y), the restored image f(x,y) and the degradation function h(x,y):J=∥g(x,y)−h(x,y)*f(x,y)∥2 
(* indicates a convolution integral). The above expression means that the evaluation amount J is given by the size of a difference between the image h(x,y)*f(x,y) obtained by causing the degradation function h(x,y) to operate on the restored image f(x,y), and the actual blurred image g(x,y). If the restored image has been correctly restored, h(x,y)*f(x,y)=g(x,y) and the evaluation amount is 0, in theory. It is assumed that the lower the evaluation amount J is, the better the restored image f(x,y) has been restored. In the steepest descent method, the iterative computation is repeated until the size of VJ which is a gradient of this evaluation amount J, that is, the square of norm of VJ becomes smaller than or equal to the threshold, and when it becomes smaller than or equal to the threshold, the iterative computation is finished and the restored image f(x,y) is obtained. It should be noted that although the case of the general inverse filter has been illustrated for the evaluation amount J, it is of course applicable to arbitrary evaluation amounts for evaluating the difference between g(x,y) and h(x,y)*f(x,y), also including the evaluation amounts in the Tiknohov-Miller normalization method and the Bayes restoration method for calculating the difference between g(x,y) and h(x,y)*f(x,y).
Now returning to FIG. 18, the image restoration processing apparatus determines whether or not t exceeds the threshold Thr (S26). If t exceeds the threshold Thr, since the restoration is regarded as still insufficient, VJ is multiplied by the convergence parameter ε (S28), and a new restored image is generated by subtracting ε VJ from the restored image (S30). The process of step S18 to step S30 is iterated until t becomes less than or equal to the threshold Thr.
If t becomes less than or equal to the threshold Thr (the determination result at step S26 is positive “Y”), or even if t does not become less than or equal to the threshold Thr, if the maximum iteration count has been reached (the determination result at step S18 is negative “N”), the process is terminated.
The above described method is the method of performing the image restoration by using the iterative computation with the steepest descent method.
It should be noted that the principle of the image restoration process involving the iterative computation is described, for example, in Non-patent Document 1 (see M. Elad and A. Feuer; Super-Resolution of An Image Sequence—Adaptive Filtering Approach; Technion—Israel Institute of Technology, 4 Apr. 1997).